Bayesian Optimization Using Sequential Monte Carlo

We consider the problem of optimizing a real-valued continuous function f using a Bayesian approach, where the evaluations of f are chosen sequentially by combining prior information about f , which is described by a random process model, and past evaluation results. The main difficulty with this approach is to be able to compute the posterior distributions of quantities of interest which are used to choose evaluation points. In this article, we decide to use a Sequential Monte Carlo (SMC) approach. 1 Overview of the contribution proposed We consider the problem of finding the global maxima of a function f : X → R, where X ⊂ R is assumed bounded, using the expected improvement (EI) criterion [1, 3]. Many examples in the literature show that the EI algorithm is particularly interesting for dealing with the optimization of functions which are expensive to evaluate, as is often the case in design and analysis of computer experiments [2]. However, going from the general framework expressed in [1] to an actual computer implementation is a difficult issue. The main idea of an EI-based algorithm is a Bayesian one: f is viewed as a sample path of a random process ξ defined on R. For the sake of tractability, it is generally assumed that ξ has a Gaussian process distribution conditionally to a parameter θ ∈ Θ ⊆ R, which tunes the mean and covariance functions of the process. Then, given a prior distribution π0 on θ and some initial evaluation results ξ(X1), . . . , ξ(Xn0) at X1, . . . , Xn0 , an (idealized) EI algorithm constructs a sequence of evaluations points Xn0+1, Xn0+2, . . . such that, for each n ≥ n0, Xn+1 = argmax x∈X ρ̄n := ∫ θ∈Θ ρn(x; θ)dπn(θ) , (1) where πn stands for the posterior distribution of θ, conditional on the σ-algebra Fn generated by X1, ξ(X1), . . . , Xn, ξ(Xn), and ρn(x; θ) := En,θ((ξ(Xn+1)−Mn)+ | Xn+1 = x) is the EI at x given θ, with Mn = ξ(X0) ∨ · · · ∨ ξ(Xn) and En,θ the conditional expectation given Fn and θ. In practice, the computation of ρn is easily carried out (see [3]) but the answers to the following two questions will probably have a direct impact on the performance and applicability of a particular implementation: a) How to deal with the integral in ρ̄n? b) How to deal with the maximization of ρ̄n at each step? We can safely say that most implementations—including the popular EGO algorithm [3]—deal with the first issue by using an empirical Bayes (or plug-in) approach, which consists in approximating πn by a Dirac mass at the maximum likelihood estimate of θ. A plug-in approach using maximum a posteriori estimation has been used in [6]; fully Bayesian methods are more difficult to implement (see [4] and references therein). Regarding the optimization of ρ̄n at each step, several strategies have been proposed (see, e.g., [3, 5, 7, 10]). This article addresses both questions simultaneously, using a sequential Monte Carlo (SMC) approach [8, 9] and taking particular care to control the numerical complexity of the algorithm. The main ideas are the following. First, as in [5], a weighted sample Tn = {(θn,i, wn,i) ∈ Θ × R, 1 ≤ i ≤ I} from πn is used to approximate ρ̄n; that is, ∑I i=1 wn,i ρn(x; θn,i) →I ρ̄n(x). Besides, at each step n, we attach to each θn,i a (small) population of candidate evaluation points {xn,i,j , 1 ≤ j ≤ J} which is expected to cover promising regions for that particular value of θ and such that maxi,j ρ̄n (xn,i,j) ≈ maxx ρ̄n(x). 2 Algorithm and results At each step n ≥ n0 of the algorithm, our objective is to construct a set of weighted particles Gn = { ( γn,i,j , w ′ n,i,j ) , γn,i,j = (θn,i, xn,i,j) ∈ Θ ×X, w ′ n,i,j ∈ R , 1 ≤ i ≤ I, 1 ≤ j ≤ J }